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ABSTRACT 

We investigate the dynamical effect of the turbulence in baryonic intergalactic medium 
(IGM) on the baryon fraction distribution. In the fully developed nonlinear regime, the IGM 
will evolve into the state of turbulence, containing strong and curved shocks, vorticity and 
complex structures. Turbulence would lead to the density and velocity fields of the IGM to be 
different from those of underlying collisionless dark matter Consequently, the baryon frac- 
tion fb will deviate from its cosmic mean jcosmic study these phenomena with simulation 
samples produced by the weighted essentially non-oscillatory (WENO) hybrid cosmological 
hydrodynamic/N-body code, which is effective of capturing shocks and complex structures. 
We find that the distribution of baryon fraction is highly nonuniform on scales from hundreds 
kpc to a few of Mpc, and varies from as low as 1 % to a few times of the cosmic mean. We 
further show that the turbulence pressure in the IGM is weakly scale-dependent and compa- 
rable to the gravitational energy density of halos with mass around 10^^ h~^ Mq. The baryon 
fraction in halos with mass equal to or smaller than 10^^ h~^ Mq should be substantially 



lower than /, 



cosmic 
b 



Numerical results show that is decreasing from 0.8/; 



cosmic Jj^Jq mass 



scales around 10^^ Mq to 0.3/™*^™"^ at 10^^ h^^ Mq and shows further decrease when 
halo mass is less than 10^^ h~^ Mq. The strong mass dependence of /{, is similar to the ob- 
served results. Although the simulated /{, in halos are higher than the observed value by a 
factor of 2, the turbulence of the IGM should be an important dynamical reason leading to the 
remarkable missing of baryonic matter in halos with mass < 10^^ h^^ Mq. 
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1 INTRODUCTION 

In the concordance ACDM universe, the baryonic gas (IGM) traces 
the collisionless cold dark matter in the linear regime of gravita- 
tional clustering. However, observation shows that these two com- 
ponents do not related to each other by a linear mapping on scales 
up to a few Mpc. Most attention on the reason of the IGM-cold 
dark matter separation on such scales has been drawn to the ther- 
mal property of baryonic gas, mainly radiative heating and cooling, 
such as photo-heating, feedbacks from star formation and accretion 
to black holes (e.g., Valageas & Silk 1999; Tozzi & Norman 2001; 
Voit et al. 2002; Zhang & Pen 2003; Xue & Wu 2003). 

The hydrodynamical origin of the IGM-dark matter separa- 
tion has been noticed in the early studies of large scale struc- 
ture foiTnation (Shandarin & Zel'dovich 1989). Successive studies 
found that in the weak and moderate nonlinear regime of cluster- 
ing, the evolution of collisionless particle can be described by the 
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Zel'dovich approximation, while the dynamics of baryonic compo- 
nent is sketched by the adhesion approximation with dissipation of 
kinetic energy, or by a random force driving Burgers equation (Gur- 
batov, Saichev, & Shandarin 1989; Berera & Fang 1994; Vergas- 
sola et al. 1994; Jones 1999; Matarrese & Mohayaee 2002; Pando, 
Feng, & Fang 2004). When the Reynolds number is large. Burgers 
fluid will be turbulent, consisting of shocks (Lassig 2000). Burgers 
shocks occur not only in high, but also in middle and even low den- 
sity regions. These shocks could contribute partly to the deviation 
of the velocity and density fields of the IGM from that of cold dark 
matter (Pando et al 2004; Kim et al 2005). 

A new progress shows that in the highly developed nonlin- 
ear regime the velocity field of the IGM is no longer potential, but 
dominated by vorticity on scales from one and a half hundred kpc 
to a couple of Mpc (Zhu, Feng & Fang 2010). Oblique shocks will 
foiTn in inhomogeneous baryonic fluid and act as the source of vor- 
ticity of the IGM. These results give supports to the scenario that in 
the nonlinear regime the baryonic fluid is in the state of fully devel- 
oped turbulence, which can be characterized by the She-Leveque' 
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(SL) scaling (He et al. 2006) and the log-Poisson hierarchy (Liu 
& Fang 2008). Although tiny vorticity can also be generated in 
dark matter field during shell crossing, the velocity field of dark 
matter is still dominated by potential motion on scales around 1 
Mpc (Pichon & Bcrnardcau 1999; Pucblas & Scoccimarro 2008). 
Consequently, the velocity and density fields of the IGM would de- 
part from those of underlying coUisionless dark matter on scales 
below a few Mpc. Some features predicted from the fully devel- 
oped turbulence have been found to be consistent with observation, 
such as the log-Poisson non-Gaussianity of Lya transmitted flux 
of quasar's absorption spectrum (Lu , Chu, & Fang 2009; Lu et al. 
2010), the scaling relations among the X-ray luminosity, tempera- 
ture and SZ effect of clusters (Zhang et al 2006; Yuan et al. 2009) 
and the turbulence broadening of HI and Hell Lya absorption lines 
(Zheng et al. 2004; Liu et al. 2006). 

In this paper, we extend these studies to the baryon fraction 
fb, which is defined as the mass ratio between the baryonic and 
total matter in a system. In the concordance ACDM universe, the 
cosmic mean of baryon fraction is = Q.i,/Q.m ~ 0.17 ± 

0. 01, where Qb and flm are the mean mass density parameters of 
the baryonic matter and total matter respectively (Dunkley et al. 
2009; Komatsu et al. 2009). The baryon fraction in gravitational 
bound objects is found to be lower than the cosmic mean, known 
as the missing baryon problem. For galaxy clusters and groups, the 
baryon fraction is smaller than the cosmic mean by a factor of 2 - 4 
(Ettori 2003; Giodini et al. 2009; Dai et al. 2010), while can be as 
low as about 1% of the cosmic mean for dwarf galaxies(McGaugh 
et al. 2010). Generally, the baryon fraction decreases monotonically 
with decreasing mass of collapsed halos. This problem has also 
been seen in the study of galaxies abundance. The ACDM model 
predicts too many low mass dark matter halos in comparison with 
estimation from the observed luminosity function of galaxies. It 
implies that the star formation rate in low mass dark matter halos is 
substantially low. The baryon content residing in virialized objects 
with small mass should be much less than that given by the cosmic 
mean. Additional physics are needed to keep the baryonic matter 
from overcooling (White & Frenk 1991). 

Many kinds of mechanisms have been proposed to prevent 
most of the baryon falling into dark matter halos. Feedback from 
supemovae and AGN and photo-heating are well investigated. 
Feedback from massive stars and SNe were once believed to be 
able to driving out the collapsed gaseous baryon and hence suppress 
the starformation(Dekel & Silk 1986; White & Frenk 1991; Kauff- 
mann et al. 1999). Theoretically, simulations and observations have 
shown that this mechanism is unlike to work while taking appro- 
priate feedback efficiency, mass loss rate(e.g. Mac Low & Ferrara 
1999; Benson et al. 2003), except for low mass halos. 

In contrast, mechanisms in which the missing baryon are in- 
hibited to fall into the protogalactic halos in the very beginning, 

1. e., never fell into, are much preferred (Mo et al. 2005; Anderson 
& Bregman 2010). Photo-heating fulfills this category which adopt 
a photo-ionizing field to reheat the baryon around progalactic halos 
and hence keep them from collapse into potential wells(e.g. Efs- 
tathiou 1992; Gnedin 2000; Benson 2002). The efficiency of this 
model, however, was later found to be largely limited, only works 
for halos below 10^" M0 (Hoeft et al. 2006), because the central 
self-shielding will delay the photo-evaporation. 

In our view, IGM turbulence should be an important reason of 
the missing baryon problem. The effects of turbulence on the prop- 
erties of clusters have been extensively studied by numerical simu- 
lations and theoretical works(e.g. Norman & Bryan 1999; Dolag et 
al. 2005; Vazza et al. 2009, 2010; Bums, Skillman& O'Shea 2010; 



Ruszkowski & Oh 2010). Yet, we will focus on the dynamical and 
statistical effect of turbulence on the spatial distribution of baryon 
fraction. That is, the turbulence of IGM is treated as a basic envi- 
ronment factor of gravitational clustering. The baryon missing in 
virialized objects should be a result of the formation and evolution 
of inhomogeneity of spatial distribution. One can then explain 
the baryon missing by considering the inhibition of gravitational 
collapse by turbulence pressure (Chandrasekhar 1951a, 1951b). 

This paper is organized as follows. §2 addresses the theoretic 
background of turbulence in the IGM and its effect on the gravita- 
tional collapsing of baryonic matter. §3 gives the simulation method 
of producing samples. In §4, we study the properties of the non- 
uniform distribution of the baryon fraction, and its relation with 
turbulence. §5 presents the baryon fraction in collapsed halos. We 
discuss our results and compare them with previous numerical stud- 
ies in §6. Finally, the conclusions are given in §7. 



2 THEORETICAL BACKGROUND 
2.1 Turbulent IGM 

Hierarchical structure formation process will introduce shocks(Ryu 
et al. 2003; Pfrommer et al. 2006; Vazza et al. 2009). The shock 
wakes and shear flows, arose from gas accreting into pancakes, fil- 
aments and halos, protogalactic and collapsed objects moving in 
complex structures and gaseous structure colliding and merging, 
will interplay with each other and drive instability, result in the on- 
set of turbulence. Fully developed turbulence consists of eddies on 
various scales and will cascade the kinetic energy of largest eddies 
down to smaller one (Lin 1966; Shu 1992). 

More specifically, the evolution of baryonic fluid in the mod- 
erate nonlinear regime can be characterized by random force driv- 
ing Burgers equation (Gurbatov et al. 1989; Berera & Fang, 1994; 
Jones, 1999; Matarrese & Mohayaee, 2002). When the effective 
Reynolds nvunber is large, burgers fluid will be turbulent, consist- 
ing of shocks and complex structures on various scales(Polyakov, 
1995; Lassig 2000; Boldyrev et al. 2004). As dark matter is not 
influenced by Burgers turbulence, the IGM velocity field will dy- 
namically decouple from the dark matter field on scales larger than 
the Jeans length once Burgers turbulence is developed. This will 
lead to the deviation of fb from in low and high density 

areas. 

The velocity field of baryonic fluid keeps irrotational in the 

moderate nonlinear regime. Latter on, the velocity field of IGM 
will no longer be potential dominated, because the vorticity can 
be effectively generated by oblique shocks due to baroclinic insta- 
bility (Zhu et al 2010). When shocks propagate in inhomogeneous 
medium, they will generally evolve into oblique or curved shock 
waves (e.g. Landau & Lifshitz 1987), acting as the source of vor- 
ticity (e.g. Picone et al. 1984; Emanuel 2000). Once triggered, the 
vorticity can be self-amplified by the nonlinearity of hydrodynam- 
ics (see §2.3). Finally the IGM evolves into a fully developed turbu- 
lent state which satisfies the She-Leveque scaling and log-Poisson 
hierarchy (He et al. 2006; Liu & Fang 2008; Lu et al 2009, 2010). 
More details about the development of turbulence in the IGM along 
cosmit evolution can be found in (Zhu et al. 2010; Fang & Zhu 
2011). 
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2.2 Turbulence pressure and vorticity of IGM 

The Jeans length of self-gravitational clustering in a statistical uni- 
form gas with density p and speed of sound Cs with or without 
turbulent gas motion has been studied by Chandrasekhar (1951a, 
1951b). In the absence of turbulence, the Jeans length is Aj = 
Cs yvr/Gp, i.e. gravitational clustering can occur only if the per- 
turbation scale is larger than Aj. The effect of turbulent motions on 
the clustering is to replace Cs by an effective sound speed 
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where {v'^) measures the velocity fluctuations of turbulence. The 
Jeans length is increased by the random velocity field of turbu- 
lence, which plays the similar role as the randomly thermal motion. 
The random motion of the turbulence provides an extra pressure, 
Ptub = (pf called turbulence pressure. The turbulence pressure, 
joining the thermal pressure, will slow down and even halt the IGM 
falling into a gravitational well. 

Considering the gravitational collapsing on scale R will not be 
affected by the velocity dispersion on scales that larger than R, the 
velocity fluctuations with wave-numbers k < 2tt/R will not resist 
gravitational collapsing on scales larger than R. Consequently, the 
effect of turbulence pressure on gravitational collapsing on scale R 
could be estimated by (Bonazzola et al 1987) 



/•Kmax 

PtuT = / E(k)dk, 

J max(fcp{ ,fcmin) 



(2) 



where E{k) is the power spectrum of kinetic energy density 
(l/2)p(r)v^(r) of the turbulence. The upper limit of the integral 
eq.(2) fcmax = 27r//diss corresponds to the minimal scale /diss be- 
low which the turbulence is dissipated. The lower limit of the inte- 
gral eq.(2) is the maximum of fen, and fcmin, where fca = 2ti / R and 
femin matches the upper scale of turbulence. Obviously, the turbu- 
lent pressure ptur is dynamical, not thermal. It can be comparable 
to the thermal pressure of the IGM, especially in regions the tem- 
perature of IGM is not very high. 

When we estimate the turbulence pressure with eq.(l) or (2), 
we should separate the velocity field v(r, t) into two components: 
one is the random motion of turbulence and the other is the bulk 
velocity. The latter mainly depends on the gravitational potential, 
and does not contribute to the turbulence pressure. In some algo- 
rithms, the bulk velocity is identified as the mean velocity within 
a box, while the fluctuation with respect to the mean velocity is 
supposed to be the turbulent motion. As the box size is selected by 
hand, these algorithms might contain non-ignorable system bias. A 
better way to pick up the random motion of turbulence bases on the 
vorticity of velocity field. 

The vorticity field uj(Y, t) of the velocity field v of the IGM 
is defined by Ui = {l/2){diVj — djVi), i.e., a; = V x v, where 
i = 1, 2, 3. The dynamical equation of the vorticity w is (Zhu et al. 
2010) 



Doj ^ ^ 1 „^ 

— — = dtu + -V • Vw 
Dt a 

= -(S ■ oj — dcj + X Vp — duj), 
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(3) 



where p is the pressure of the IGM, d = diVi is the divergence 
of the velocity field and a{t) is the cosmic factor. Tensor S, de- 
fined as Sij = {l/2){diVj + djVi), is called strain rate, and 
[S • w]i = SijOjj. A remarkable property of eq.(3) is the absence of 
the gravity term . This point is expected, as gravity is curl-free in 



nature. In other words, vorticity is fully given by the nonlinearity 
of hydrodynamics and hence is effective to measure the turbulence. 

An important property of the vorticity is that for a fully de- 
veloped turbulence the power spectra of the vorticity field Pu(k) 
and the velocity field Pv(k) should satisfy Puj{k) — k^Pc{r) 
(Batchelor 1959, 2000). This relation can be used to determine the 
wavenumber fcmin , i.e. the upper scale of fuUy developed turbu- 
lence. We use this criterion to ascertain fcmin to calculate the inte- 
gral eq.(2). 



2.3 Dynamical effect of turbulence on IGM gravitational 
clustering 

The dynamical effect of turbulence on the gravitational clustering 
of the IGM can be seen with the time evolution of the irrotational 
component of velocity field, i.e. the divergence of velocity field 
d = diVi, which follows (Zhu et al 2010) 
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where ptot is the total mass density including both cold dark and 
baryonic matter with po being its mean value. As negative diver- 
gence means the increase of density, positive AnG[ptot — po)/o, 
will lead to clustering, while negative one result in anti-clustering. 
The right hand side of eq.(4) can be used to compare the effects of 
hydrodynamical terms on clustering with that of gravity. 

We first identify the physical meaning of the hydrodynamical 
terms on the right hand sight of eq.(4) by considering incompress- 
ible fluid, i.e. assuming ptot to be a constant ptot = po- In this 
case, the term of gravity (47r/a) (ptot — po) = 0, and Vp = 0. The 
continuity equation yields d = Q, and Eq.(4) gives 



(5) 



V p = -p ( Sij Sij - -uj' 



This is a typical Poisson equation for the scalar field p. Analogous 
to the field equations in electrostatics, the term p[SijSij — 1/2lj^] 
on the right hand side of eq.(5) plays the role of the "charge" of 
a pressure field. Positive "charge" produces an attraction force, 
while negative "charge" yields a repulsive force. Back to eq.(4), 
p[SijSij — l/2t<;^] also plays the role of nonthermal pressure 
of turbulence. In regions with negative "charge", i.e. [Sij Sij — 
(1/2)0)^] < 0, the turbulent fluid will prevent the gravitational 
clustering. 

The sign of the "charge" is actually determined by the levels 
that turbulence has developed. From the definition of vorticity and 
strain rate, we have 

1 2 1 

-UJ - S^jS^j = -[{diVj){diVj) - ?,{djV{){diVj)]. (6) 

For a Gaussian velocity field, {3{djVi){diVj)) = {{diVj){diVj)), 
the net effect of velocity fluctuations on the IGM collapsing is 
null in average. However, for a non-Gaussian velocity field, it 
can be either positive or negative, depending on the property of 
the velocity field. For a homogeneous and isotropic turbulence , 
({diVj){djVi)) = (Batchelor 1959), the sign of the "charge", 
p[SijSij — l/2a;^], is always negative. A fully developed turbulent 
flow will generally prevent the IGM from gravitational collapsing. 
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The term V^p/p of eq.(4) relates to the hydrostatic pres- 
sure while does not consider fluid compressibility, which has been 
presented in eq.(5). It is mostly negative in overdense collaps- 
ing regions, and will resist upon gravitational collapse. The term 
— (Vp) ■ {yp)/f? on the right hand side of cq.(4) result from the 
compressibility of the IGM. Its value would be negative when the 
density-pressure relation is a power law p oc and 7 > 0, and 
then also plays the role of resisting gravitational collapsing. 

The physical meaning of the term (Vp) • (Vj3)/p^ can be 
shown by the ratio between (Vp) ■ (Vp)/ and the gravity term 
47rG(ptot - Po). The ratio is roughly equal to ~ (tinfau/^sound)^, 
where tinfaii ~ {GpY^^ is the free falling time scale, tsound ~ 
R/cs, Cs ~ (Vp/Vp)^^^, is the dynamical time scale on a spatial 
scale R collapsing. When (iinfaii/isound)^ is larger than 1.0, the 
free falling time scale is larger than the dynamical time scale, and 
in result, the collapsing on scale R is significantly prevented. 



3 METHOD 

As mentioned above, the turbulent IGM contains curved shocks, 
vortices, and other discontinuities. To study the effect of turbulent 
IGM with simulation, the algorithm should be qualify for captur- 
ing these complex structures. We take advantage of the develop- 
ment of the Weighted Essentially Non-oscillatory (WENO) method 
(Shu 1999). The WENO schemes have been widely used in various 
fields, such as high Reynolds number compressible flows (Zhang et 
al. 2003), high Mach number jets (Carrillo et al. 2003), magneto- 
hydrodynamics (Jiang & Wu, 1999) and hypersonic boundary layer 
(Rehman et al 2009). The shock capturing algorithm with WENO 
scheme has past many test, including shock-boundary layer in- 
teraction (Lagha et al. 2009), shocks in high-speed flows (Mar- 
tin, Piomelli & Candler 2000), shock vortex interaction (Grasso & 
Pirzzoli 2000a, 2000b) and shock-turbulence interaction (PirozzoU 
2002). An updated review of the WENO method is given in Shu 
2009. 

In the context of cosmological hydrodynamical simulation, 
the WIGEON code is based on Eulerian description of hydrody- 
namics with 5th order WENO finite difference scheme and particle 
mesh (PM) method for dark matter particles (Feng, Shu, & Zhang 
2004). The WIGEON code can reproduce commonly accepted re- 
sults such as the relation between IGM temperature and density and 
the component of WHIM (He, Feng, & Fang 2004). In the same 
time, it also reveals a series of turbulent behavior of the IGM (He 
at al 2006; Liu & Fang 2008; Lu et al. 2009, 2010). These features 
of this code fit well with our goal. 

The cosmological parameters are taken to be 
(a^,aA,ft,(T8,fib,ns,«r-e) = (0.274, 0.726, 0.705, 0.812, 
0.0456, 0.96, 11.0) (Komatsu et al, 2009). The simulation run in a 
periodic cubic box of size 25 Mpc since redshift z = 99 with 
a 512* grid and an equal number of dark matter particles, giving 
a mass resolution of 1.04 x 1O^M0. A uniform UV background 
of ionizing photons is added at Zre to mimic the reionization. 
Radiative cooling and heating are followed as Theuns et al. 
(1998) with the primordial composition X = 0.76, Y = 0.24. 
Star formation and its feedback are not added in our simulation, 
because the resolution we used makes it hard to implant these 
processes appropriately. Snapshots are outputted at redshifts 
z = 11.0, 6.0, 4.0, 3.0, 2.0, 1.0, 0.5, 0.0. To study the convergence 
of numerical results, we also run a simulation with 256^ grid and 
an equal number of dark matter particles. 

The velocity field of simulation samples has been used to show 



the development of vorticity and turbulence in Zhu et al.(2010). To 
locate shocks post simulation we use a algorithm also based on 
WENO kernel, combining with conditions from Ryu et al.(2003) 
(see Appendix). To construct the density field of dark matter on 
grids, we assign the mass of dark matter particles onto grids using 
the Triangle-Shaped-Cloud (TSC) method. In order to minimize 
the system bias of assignment in calculating the baryon fraction on 
grids, we smooth over the density fields of baryonic and dark matter 
on grids separately using the same smooth window with radius of 
one grid. 

To identify halos of dark matter and their radius r2oo, we use 
the same process as Grain et al. (2007). Halos are identified by 
two methods: a.) friends-of-friends (EOF) method with a linking 
length parameter 0.2; b.) HOP method (Eisenstein & Hut 1998) 
with default parameters for group searching and (5outer = 80, 
^saddle = 200 and 5peak = 240 for group merging. We study 
only halos consisting of no less than 2000 dark matter particles at 
z = 0. For each halo, we can then find the center and radius of 
a sphere, in which the mean mass density is equal to 200pcrit(2), 
where pcrit(^;) = 3H^{z)/8itG is the critical density at redshift z. 



4 SPATIAL DISTRIBUTION OF BARYON FRACTION 
4.1 A Slice 

Figure 1 presents an example of the spatial distributions of bary- 
onic and dark matter density in a slice of 25x25x0.2 fr'-^ Mpc^ at 
redshift z = 0. The density fields of dark matter and baryonic fluid 
of Figure 1 display the typical sheets-filaments-knot structures on 
the cosmic scales. 

The spatial distributions of shock, temperature, vorticity and 
normalized baryon fraction F[, = /[,/ on the same slice are 
shown in Figure 2. The shocks do not always follow the filament 
and sheet structures of the matter density fields. Shocks can also be 
formed in low density areas, as demonstrated in other simulations 
(e.g. Ryu et al. 2003; Kang et al. 2007). 

The vorticity field in Figure 2 is described by the dimension- 
less quantity ojt, where t is the cosmic time, ojt accounts the num- 
ber of rotation of the vorticity w in the age of universe. The distri- 
bution of vorticity also doesn't always follow the filament and sheet 
structures of the matter density fields, but has cloud-like structures 
(Zhu et al 2010). 

Although curved shocks are the sources of vorticity, the spatial 
distribution of vorticity field does not show the same configuration 
as the shocks. It is because the spatial transfer of vorticity is given 
by the velocity field [eq.(3)], while the case with shock front is dif- 
ferent. Once the vorticity appears, it will depart from their source- 
curved shocks and spreads over the space along with the velocity 
field (Batchelor 2000). 

In the plot of vorticity distribution, we add FOE identified ha- 
los, which are marked as solid circles. The radius of circles are 
enlarged to 5 times of the real halos radius. A remarkable feature is 
that all the halos arc located in the area of Igi^t > 1. It means that 
all the gravitational bounded halos are formed in the environment 
of turbulent IGM. As discussed in §2.2, the gravitational clustering 
should be affected by the turbulence pressure. 

The spatial distribution of the normalized baryon fraction Fb 
is highly non-uniform. The value of Ft spreads over three orders of 
magnitude from Ig ~ — 2 to 1, indicating that the separation 
of baryonic matter from dark matter is significant developed dur- 
ing the nonlinear evolution. The deviation of baryon fraction from 
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Figure 1. Distributions of dark matter density (left); baryon density (right) in a slice of 25x25x0.2 h ^ Mpc^ at redshift z = 0. 



the cosmic mean is not only occurred in or around massive halos, 
but also in the low and moderate density areas. The spatial dis- 
tribution of baryon fraction does not display the sheets-filaments- 
knot structures as the density fields of dark and baryonic matter do. 
Nevertheless, the area of high Pcdm are surrounded by high baryon 
fraction area (dark blue), indicating that the processes of separat- 
ing baryonic matter from dark matter would be potent around high 
density areas. The sizes of some region in which IgFt < —0.5, 
corresponding to Ft < 0.3, is much larger than the resolution, 
demonstrating that the existence of regions with Ff, < 0.3 is ro- 
bust. 



4.2 Comparing the effects of turbulence and gravity on 
clustering 

As mentioned in §2.3, the term (47rG/a)(ptot — po) on the right 
hand side of eq.(4) measures the gravity effect leading to clus- 
tering (ptot > po), or anti-clustering (ptot < po)- The first 
three terms on the right hand side of eq.(4) describe the gas dy- 
namical and thermal effects, namely, the vorticity and strain rate 
—uj'^/2 + SijSij, the ordinary pressure V'^p/p and multiphase- 
related term — (Vp) ■ (Vp) / p^. We now use simulation sample to 
study the general properties of these terms. 

Figure 3 presents a comparison between ( Vp) ■ ( Vp) /p^ and 
47rG(ptot — Po) in cells randomly selected from the simulation 
samples at 2 = 0. We chose only the cells of 47rG(ptot — po) > 0, 
which drives collapsing. The variables used in Fig. 3, (l/p^)(Vp) ■ 
(Vp)t^ and 47rG(p tot — po)i^, are dimensionless, where t is the 
age of the universe. 

We can see from Figure 3, when 47rG(ptot — po) < 50, 
more than a half of the data points having (Vp) • (Vp)/p'^ > 
47rG(ptot — po) and few points have (Vp) • (Vp)/p^ < 0, which 
corresponds to the so-called inverse density-pressure(temperature) 
relation. Therefore, the term (Vp) • (Vp)/p^ enhanced by shocks 
and complex structures can effectively resist gravitational col- 
lapsing till 47rG(ptot — Po) ~ 50. For deep gravity wells with 
47rG(ptot — po) > 50, most data points show that gravity gener- 
ally is larger than (Vp) • iWp)/ p^ . Namely, only in deep gravita- 
tional wells the motion of baryonic matter will be dominated by the 
gravity. 

Figure 4 compares — SijSij vs. (Vp) • iWp)/ p^ (left 

panel), and (l/2)aj'^ — SijSij vs. (l/p)V'^p (right panel) at ran- 
domly selected cells. We use the dimensionless variables by multi- 



plying with . The left panel of Figure 4 shows that the magnitude 
of a;^/2 — SijSij is comparable with (Vp) ■ iyp)/ p^. While in 
the right panel, it exhibits that term (l/2)a;^ — SijSij is stronger 
than the term (1/p) V'^p. On average, the term (1/p) V^p is smaller 
than (Vp) ■ (Vp)/p^ 

In a word, in terms of prevention of gravitational collaps- 
ing, the effect of vorticity and strain rate is comparable with 
(Vp) • (Vp)/p^ and dominates over (l/p)V^p. Both of the latter 
two terms are thermal pressure p related, one can regard vorticity 
and strain rate as an effective thermal pressure. Figure 2 shows that 
in most regions IGM temperature are equal to or less than 10^ K. 
Giving the comparison in the last paragraph, the effective thermal 
pressure of turbulence would be as strong as the thermal pressure 
with temperature ~ 10^~® K. The turbulence pressure level given 
by simplified estimation is consistent with the more delicate one 
using the energy density of turbulence in Zhu et al. (2010). 

4.3 Correlations between baryon fraction and turbulence 

As vorticity is an important dynamical factor to cause a low baryon 
fraction, one can expect a correlation between the baryon frac- 
tion and the vorticity. Figure 5 gives the probability density func- 
tion(PDF) of lg[(cjf)V2] in various ranges of the baryon fraction 
Fb. Figure 5 indeed shows that for cells have small Fb, the PDF at 
the side of \g[{u]tY /2] > 2 are higher The cells with more serious 
baryon missing have larger vorticities in average. 

The similar PDF of Ig | ( 1/p) {V^p)t^ \ is also shown in Figure 
5. The cells with more significant baron missing generally corre- 
sponds to large |(l/p)(V'^p)f'^| as well. Although the two panels of 
Figure 5 look similar, the values of {ujtf /2t^ and |(l/p)(V^p)t^i 
actually cover different magnitude range. The former can be as 
large as lO'', while the later is about 10^. This is consistent with 
the right panel of Figure 4. The vorticity term (a;t)^/2 is larger 
than the term (l/p)(VV)t^. 

Figure 6 plots the Fb vs. mass density of dark matter in ran- 
domly selected cells in which (1/2)lj^ — SijSij is positive and 
{1/2)lj^ — SijSij dominant over |47rG(ptot —po)|- It shows clearly 
that the higher the density pcdm the lower the baryon fraction Fb- 
Since high density pcdm means high |47rG(ptot — Po)|, Figure 6 
indicates that as long as the turbulence of the IGM is dominant, the 
baryon fraction is lower for higher Pcdm- 

We can also see from Figure 6 that a significant fraction of 
those data point which have Fb < 0.3, i.e. more than 70 percent of 



6 Zhu, Feng & Fang 



'■-.ho 


k 1 (ir.dlion 
















! t ' > 


O.v 0,.1 0.0 O.ti 1 
'1 f ' 








• 








■ ■ < 








• 

• 




• 



I 




0.2 0.4 0,6 O.S 1.0 



O.-l 0.6 




0,0 0,2 



Figure 2. The top left shows shocks with Mach number Ma ~ 1 — 10 (green), ---^ 10 — 100 (yellow) and > 100 (red). The top right is the temperature 
contours with green, navy, cyan and red representing 10^, 10^, 10^, 10^ K, respectively. The bottom left is the vorticity in the slice, in which the solid circles 
are halos identified by FOF method. The radius of each sohd circle is equal to 5 times of the halo radius. The bottom right is the normalized bai'on fraction 
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Figure 3. (Vp) • vs. 4TTG{ptot — po) of randomly selected cells at z = 0. The dotted line shows (Vp) ■ {Vp)/p^ = 47rG(ptot — po)- 



baryonic matter are missing, locate in the areas with mass density 
Pcmd ~ 1- The event of Ft < 0.3 baryon missing can occur even 
in the regions out of collapsed clumps, indicating that the baryon 
missing of viriahzed halos should be considered as result of the 
evolution of baryon distribution on sizes larger than virialized ha- 
los. 



5 BARYON FRACTION OF DARK MATTER HALOS 
5.1 Turbulence pressure and halo mass 

It has been shown that the relation P^{k) = k^Py{r) is well sat- 
isfied by the IGM velocity field in the scale range from 0.2 to 3 
Mpc at 2 = (Zhu et al 2010), on which the IGM is in the 
state of fully developed homogeneous and isotropic turbulence. The 
power spectrum of kinetic energy density of the turbulence is ap- 
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proximately E{k) oc fc^^'^. The turbulence pressure (or the kinetic 
energy density) is weakly scale-dependent ptur(i?) oc 7?"" '* and 
about the order of magnitude ~ 1 x 10~^® g cm~^ s~^. On the 
other hand, the viral temperature T200 of a halo at the virial radius 
r2()() is oc r?2QQ- The kinetic energy density of the virial motion at 
r2oo is then - 1 x 10"^'^(r2oo/10'^ ®). The mass of halos with 
T200 = 10^-^ K is about lO" M®. Obviously, the effect of 
turbulence pressure on the IGM collapsing would be comparable 
with and even larger than the gravity for halos with mass equal to 
and less than 10^^ h^^ Mq. The baryon fraction would signifi- 
cantly be reduced for halos with mass < 10^^ Mq. 

Meanwhile, the kinetic energy density of turbulence is small 
than that of the virial motion of halos with mass larger than times 
of 10" Mq. Thus, the halo mass dependence of the baryon 
massing would show two phases: when the halo mass is larger than 
times of 10^^ Mq, the baryon fraction ft, is not much affected 
by tmbulent IGM, and it will be substantially decreasing with halo 
mass when halo mass is lower than a few 10^^ Mq. 

5.2 Halo mass dependence of baryon fraction 

With the preparation of above sections, we now analysis the 
baryon fraction in gravitational collapsed halos. The accuracy of 
the baryon fraction in selected halos relies on the calculation of the 
mass of baryonic matter enclosed in the viral radius. To reduce the 
system error, we divide each halo sphere into sub-cubic cells with 
size equal to 1/50 of the sphere radius. We then use Cloud-In-Cell 
(CIC) interpolation to determine the baryon density at each sub- 
cell. The halo mass M200 and baryon mass are given by the sum of 
the mass of total matter and baryonic matter respectively in all the 
sub-cubic cells. This method may yield large uncertainty for ha- 
los with r2oo less than the size of a grid. To restrain the uncertainty, 
our statistical results mainly are on M200 > 10^^ M©, i.e. the virial 
radius is about equal to or larger than two grids. 

Figure 7 displays the baryon fraction as a function of halo 
mass at 2: = 0, where the halos are identified by FOF (left) and 
HOP (right) respectively. The results of the two panels are statisti- 
cally the identical. At a given halo mass, there are large scatters in 
the baryon fraction F},, the scatters are seen to be even more sig- 
nificant in the FOF algorithm than the HOP. This is consistent with 
Figures 6, which indicates a large scatter in Fb at a given mass den- 
sity. The distribution of Fi, is highly non-Gaussian and has a long 
tail on the side of small Fi,, conforming with statistics of turbulent 
fluid, of which the velocity and density fields are intermittent and 
their PDFs have long tails. 

Regardless of the scatters, the distribution of Fi, in Figure 7 
can not be described by a single power law, but shows two distinct 
phases in halo mass-dependence. In the mass range 10^^ — 10^^ 

Mq, Fb is weakly dependent on halo mass and takes a value 
of 80% to 90% of the cosmic mean, which is the first phase. When 
halo mass is less than < 10^^ Mq, Fb is quickly dropping off 
with decreasing halo mass, i.e. the second phase. At mass~ 10^^ 

Mq, Fb is around 0.3 of the cosmic mean. For halos with 
mass M2QQ = 3 X 10^" M©, the mean of Fb is even as low as 
about 0.1. Although the data points in the second phase shows a 
large scattering, this phase has a clear upper envelop. The strong 
downward in Fb with decreasing halo mass is evident even taking 
account of the scatters. The scatter looks very small on the side of 
high halo mass. 

To test the convergence and stability of the results of Fig. 7, we 
1.) calculate Fb within radius 2r2oo; 2.) analysis 256^ samples (§3). 
For the 256^ samples, we only calculate the Fb within radius 2r2oo 



to avoid large system error, as r2oo is under the grid size at the low 
mass end. The results are given in Fig. 8, in which the halos are 
identified by the HOP algorithm. It shows that Fb in 2r2oo is about 
the same as that of r2oo within the error bars for halos larger than 
times of lO" h'^ Ms- For halos with mass < lO" hr'^ Ms, Fb 
in r2oo shows a little lower than that of 2r2oo- This result seems to 
be reasonable, if noted that the baryon fraction of cluster increases 
with radii (e.g. Ettori & Fabian 1999, Wu & Xue 2000). 

For the 256* simulation. Ft in radius 2r2oo shows about 10% 
lower than its counterpart of the 512'' simulation at halo mass 
10^^-10" h-^ Mq. This difference probably comes from the un- 
certainty caused by the size of grids. The uncertainty of Fb caused 
by the resolution of 512* simulation should be around 10% at halo 
mass 10" h'^ Mq 

Figure 9. compares the baryon fraction as a function of halo 
mass M200 at redshifts z = 1, and 2. The halos in Figure 9 are 
identified by HOP method, which yields similar results as the FOF 
method. Clearly, it shows the scatter in Fb becomes more signif- 
icant at lower redshifts, data points have Fb < 0.2 at 2: = 2 are 
quite few, but significantly increases at 2: = 1. 

The redshift evolution of Fb is parallel to the IGM turbulence 
(Zhu et al 2010). The turbulence is fully developed on scales from 
0.2 hr^ Mpc up to a couple of Mpc since redshift z ^ 2. In the 
IGM density distribution shown in Figure 1, there are about 7.6% 
of volume with {1/2)lo'^ — SijSij > at redshift 2 = 0, while it 
is merely 2.6% at redshift z = 2. The effect of turbulence becomes 
stronger to slow down the IGM clustering at lower redshifts. 

5.3 Comparison with observation 

Many measurements on the baryon fraction of galaxies, groups and 
clusters have been done in recent years. Although it is widely ac- 
cepted that the baryon missing occurs in gravitational bound ob- 
jects, the result has not well quantitatively settled yet. For instance, 
some measurements claim that the baryon fraction in most massive, 
relaxed galaxy clusters is close to the cosmic mean (e.g. Vikhlinin 
et al. 2006; Allen et al. 2008), while someone else argues that the 
baryon fraction of clusters with mass ~ 10^"' Mq is less than 
the cosmic mean (Giodini et al. 2009). Nevertheless, observation 
shows that the baryon fraction of objects with mass < 10^* Mq 
is systematically decreasing with decreasing masses. For galaxy 
groups and galaxies, the baryon fraction is probably not higher than 
about 0.1 of the cosmic mean (e.g. Sun et al. 2009; Hoekstra et al. 
2005; Heymans et al. 2006; Mandelbaum et al. 2006; Gavazzi et al. 
2007). 

Figure 10 shows the mass dependence of baryon fraction Fb 
measured at rsoo (MaGaugh et al. 2010) and r2oo (Dai et al (2010). 
These two data show some difference at mass > 10^* Mq. It 
is probably due to the cool gas is underestimated with X-ray mea- 
surement. It may be also partially caused by the difference in red- 
shifts of the two sample sets. 

The result of MaGaugh et al.(2010) can not be fitted by a 
single power law, but show two-phase feature. In the mass range 
> 10^* hr^ Mq, Fb ~ 0.85, almost independent of object mass. 
In the mass range < 10^* h^^ Mq, the baryon fraction decreases 
significantly with the halo mass. Fb is only about 10% of 
at 10" Mq and further downs to 6% at halo mass of 3 x 10^° 
h-^ Mq. 

In Figure 10, we also plot the mean and variance of simulation 
baryon fraction Fb in r2oo for HOP halos at z=0. The simulated 
mass-dependence of Fb shows the similar trend as observation. In 
the mass range of ~ 2 x 10^^ — 10^* h~^ Mq, simulation results 



Figure 7. Baryon fraction Fi, as function of halo mass M200 at z = 0. The halos are identified with FOF method (left) and HOP method (right). 
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Figure 8. Baryon fraction Fi, as function of halo mass Af200- The Fi, is calculated within the radius r2oo (red), 2r2oo (blue) for 512'^ samples and within 
2^200 for the 256'' samples, respectively. The error bars are the rms of the scattered distribution. 



are higher than observed data by a factor of about 2, but within the 
error bars of observed data. The effect of IGM turbulence is not 
the only reason leading to the baryon missing. Yet it should be an 
important dynamical reason of the baryon missing, especially for 
halos with mass < 2 x lO" h^^ Mq. 



6 DISCUSSION 

6.1 Comparison with previously numerical studies 

The baryon depletion has attracted many studies with cosmologi- 
cal hydrodynamical simulation. As revealed in this paper, the miss- 
ing baryon could be caused by the IGM turbulence. In order to 
accurately study the effect of the turbulent IGM on baryon fraction 
distribution and hence the baryon mission , the algorithm of hydro- 
dynamical simulation should be effective to capture curved shocks, 
vortices, and intermittence of turbulence in the IGM. 

It has been shown that the smoothed particle hydrodynamics 
(SPH) method may not be able to handle shocks or discontinuities 
as well as grid method, because the nature of SPH is to smooth be- 
tween particles(Tasker et al. 2008). The SPH method is found to 
have strong damping of velocity fluctuations and fluid shear insta- 
bilities with respect to grid method (Agertz et al. 2007). Moreover, 
the artificial viscous force (dissipation) in SPH algorithm will re- 
duce the Reynolds number, and suppress the effect of turbulence 
(Dolag et al. 2005). These factors would be part of the reasons that 



some SPH simulations find the baryon fraction in halos to be gen- 
erally only 10% lower than the cosmic mean(e.g. Eke, Navarro 
& Frenk 1998; Frenk et al. 1999; Ettori 2006; Grain et al. 2007), 
where the lower in results from energy transfer during shock 
formation(Navarro & White 1993). 

The Eulerian grid method, including both fixed grid and those 
with adaptive mesh refinement(AMR), are shown to be effective to 
pick up the turbulent behavior of baryon fluid in and around galax- 
ies and clusters (e.g. Norman & Bryan 1999; Ryu et al 2008; Mol- 
nar et al 2009; Vazza et al. 2009; Bums, Skillman& O'Shea 2010; 
Ruszkowsky & Oh 2010). However, very few Eulerian simulation 
has been conducted on the turbulent behavior of IGM other than 
Zhu et al.(2010). OShea et al. (2005) investigated the baryon frac- 
tion in halos at z = 3 during their study on the performance com- 
parison of ENZO and SPH method without looking into tur- bu- 
lence. However, their small box, merely 3 Mpc, and relative poorly 
resolved dark matter halos, typically tens to hundreds parti- cle 
in a halo, makes a meaningful comparison between their results 
and those presented in this paper unavailable. The energy cascade 
would be suppressed in small box, which would significantly re- 
strain the development of turbulence in the IGM. Very recently, the 
turbulence of IGM is studied with grid based AMR code ENZO in- 
cluding a subgrid scale model for small scale unresolved turbulence 
in lapichino et al. (201 1), which yields many results similar to Zhu 
et al (2010). However, statistical study on the baryon fraction is not 
carried out. 
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Figure 9. Baiyon fraction at 2r2oo vs. halo mass M200 at 2: = l(Ieft) and z = 2(right). The halos are identified with HOP method. 
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Figure 10. Baryon fraction vs. mass of objects given by the data of McGaugh et al. (2010) (blue); Dai et al (2010) (grey region), and simulation result of 
at r200 of the HOP halos at z=0 with a bin of 0.25 dex in mass. 



The WENO algorithm embeded in our simulation uses a con- 
vex combination of all the stencils in the reconstruction procedure, 
where each stencil is assigned with a nonlinear weight depending 
on the local smoothness. This algorithm provides a high order accu- 
rate way to capture the non-oscillatory property near strong discon- 
tinuities as well as complex smooth solution features. This feature 
makes the WENO scheme has gained rapid popularity in the sim- 
ulation of complex-structure of hydrodynamical field (Shu 1999, 
2009). The 5-th order WENO scheme used in this work makes it 
possible to follow high Reynolds number baryonic fluid. On the 
other hand, like other high order non-artificial-viscosity scheme, 
e.g. piecewise perturbation method (PPM), with fixed grid, the 
computation time needed is much more than the SPH and AMR 
method, which would become more obvious in higher resolution 
work. 



6.2 Systematic effects 

To estimate the effect of turbulence of IGM, a key factor is the 
upper scale of fully developed turbulence. In some algorithm, the 
energy of the turbulence is calculated with the fluctuations of veloc- 
ity with respect to the mean velocity in cells with selected size. An 
underestimated cell size generally leads to undervalued turbulence 
(Dolag et al. 2005). In our algorithm the upper scale of turbulence 
is determined by the relation Pvor(fc) — k^Pv(k), by which the 
ambiguity in selecting the cell size is avoided. 



Star formation and its feedback on the IGM evolution are not 
considered in our simulation. The injection of hot gas and energy by 
supernova explosions, AGN or other sources of cosmic rays is not 
followed. These factors can be properly added, if the star formation 
history and radiative transfer is well known, which in turn need very 
high resolution. As mentioned above, the influence of injecting hot 
gas and energy by supemovae is believed to yield further decrease 
of Fi, in less massive halos, while feedback from AGN might works 
in massive central galaxies. Including all these mechanisms would 
give better consistency with observation. Meanwhile, an accurate 
and robust observation result of baryon fraction in objects from 
as massive as galaxy cluster to low mass dwarf galaxy is full of 
challenge. Large scatter in the observed quantities and significant 
system bias comes from the models applied in object structure and 
radiative transfer urges much more efforts to deal with. 

The spatial resolution of our simulation is 48.8/i^^kpc, which 
is not enough to simulate virialized halos with mass less than times 
of 2 X 10^" h^^ Mq . The lower limit of the inertial range of turbu- 
lence may be affected by the grid resolution, which is already seen 
in the vorticity power spectra given by samples with different res- 
olution in Zhu et al. (2010). On the other hand, the box size of our 
simulation is 25/i^^Mpc. It may lose the effect of long wavelength 
perturbations and hence underestimate the nonlinear evolution of 
IGM velocity field. Nevertheless, we believe that the basic dynam- 
ical picture features revealed by the current samples would be valid 
when these factors are improved. 
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7 CONCLUSION 

Since gravity is of scale free, it is generally believed that in the 
scenario of hierarchical clustering, the formation and evolution of 
halos is scale free in a large range of halo mass. However, the dy- 
namics of the system consisting of two components, dark matter 
and IGM, is very different from those of one component system. In 
the nonlinear regime, the hydrodynamical nature of the IGM leads 
to the dynamical and statistical departure of the IGM from the dark 
matter. The two component system is no longer to follow the scal- 
ing of gravitational clustering of pure dark matter. This deviation is 
a reason of the baryon missing in gravitational collapsed halos. 

The dynamical equation of vorticity does not contain terms of 
gravity. Therefore, the IGM turbulence characterized by the vortic- 
ity will yield, at least, two scales, which violate the scale free of the 
gravitational hierarchical clustering: 1.) the length scale on which 
the IGM fluid has been developed to the state of fully developed 
turbulence; 2.) the mass scale on which the turbulence pressure is 
comparable with gravity of halo considered. These scales play im- 
portant role in the evolution of baryon fraction. With cosmological 
hydrodynamic simulation, we find that at 2 = the first scale is 
about 3 h^^ Mpc, and the second one is ~ 10" h"^ Mq. With 
these results, we reach to the following conclusions: 

i. The distribution of baryon fraction is highly nonuniform on 
scales from hundreds kpc to a few of Mpc, and ft varies from as 
low as 1% to a few times of the cosmic mean. 

i. The turbulence can effectively prevent the IGM from falling 
into potential wells of dark matter halos with mass ~ 10^^ h^^ 

Mq. 

iii. The ft in dark matter halos is decreasing from 0.8/^°^™"^ 
at halo mass scales around 10^^ h'^ M© to 0.3/6^^°=™'" at 10" h"^ 
Mq due to the turbulent state of the IGM. 

The estimated turbulence pressure at 2: = correspond to a 
random motion with r.m.s velocity of about 50 — 100 km in 
the scale range from hundreds of kpc and up to ~ 2 Mpc. The 
turbulent pressure is dynamical and non-thermal. When the turbu- 
lence dissipated, its kinetic energy becomes the thermal energy. It 
yields the entropy in halos (He, Feng, & Fang 2005). Therefore, 
the dissipation of turbulence actually is a mechanism of heating, 
which gives a compensation to the cooling of gas in halos. This re- 
sult is consistent with the Burgers' shock heating (He et al 2004). 
In summmary, the dynamics of turbulence can effectively affect the 
baiyon fraction of halos . 
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APPENDIX A: SHOCK CAPTURING ALGORITHM OF 

WENO SCHEME 

In the WENO scheme, the smoothness of a hydrodynamical quan- 
tities U = (p,p,v, E) in any dimension is meastired by 



(Al) 



where the definition of Ui is the same as that given in the Appendix 
ofZhuetal(2010). 

With ?7i, we can construct the one-dimensional discontinuous 
detector by 



DDFi = 



Ui+i + 2Ui + Ui-i 



(A2) 



where n = 1, 2, ... We use U = p and n = 2 to detect the discon- 
tinuous surface (Visbal & Gaitonde 2005; Lo et al. 2007). For each 
cell, we calculate DDFi in three directions. The value of DDF 
is normalized by the max of DDF among all the cells. Then, a 
threshold parameter of 1.0~* of DDFi is used to flag discontinu- 
ous zone with DDFi > 10~*. Finally, we pick up the shocks with 
Mach number Ma > 1.5 from the flagged cells by checking the 
conditions as that given by Ryu et al. 2003. 



